Exact, numerical, and mean field behavior of a dimerizing lattice in one dimension. 
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The thermodynamics and dynamics of a one dimensional dimer-forming anharmonic model is stud- 
ied in the classical limit. This model mimics the behavior of materials with a Peierls instability. 
Specific heat, correlation length, and order parameter are calculated three ways: (a) by mean field 
approximation, (b) by numerical molecular dynamics simulation, and (c) by an exact transfer matrix 
method. The single-particle spectrum (velocity-velocity correlation function) is found numerically 
and in mean field theory. 
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I. INTRODUCTION 

Polyacetylene [Q and blue bronze Q are examples of 
quasi-one-dimensional materials with broken-symmetry 
ground states pL driven by electron-phonon interaction, 
while CuGeC>3 |4j is driven by a spin-phonon interaction. 
In strictly one-dimensional (ID) models, phonon fluctu- 
ations destroy long range order at temperature T > 0. 
The electron-phonon problem has been extensively stud- 
ied and numerical studies converging towards exact an- 
swers have been made Hfjj . Most studies of the electron- 
phonon problem, especially studies of higher dimensional 
models, omit lattice dynamical fluctuations and focus on 
electronic fluctuation effects, keeping the lattice frozen. 
The present paper analyses a model with only lattice- 
dynamical fluctuations. Various authors pointed out 
that both electronic and lattice (zero-point and thermal) 
fluctuations should be taken into account to describe the 
thermodynamics of these materials. Lattice fluctuations 
work as an effective disorder on the electronic system, 
becoming stronger as the T increases jl0| . 

We take the following Hamiltonian: 




(1) 

The first two terms are a single-site (Einstein oscilla- 
tor) double-well potential with minima at ±uq where 
wo = \J ki/k,2- The next term is a second-neighbor 
spring which prefers displacements to alternate, ui = 
(— IYuq. After introducing dimensionless atom displace- 
ments x n = x n /uo, time t = t[Ki/M) 1 , and energy 
E = E/k\Uq, only one free parameter £ = k^/ki is left 
in the problem. This model was frequently studied in 
the past plHlJ]. These studies have been reviewed by 
Dieterich |15[. McKenzie considered a continuous 
version of the model. 



Here we solve this problem numerically by molecular 
dynamics (MD) simulation and compare with a mean 
field approximation (MFA) as well as an exact transfer 
matrix (TM) method. The MD simulation enables us to 
evaluate a dynamical correlation function not available 
by the TM method; we compare it with the MFA. 

II. MEAN FIELD SOLUTION. 

Exact thermodynamic calculations of Hamiltonians 
like (|l|) are possible only for a ID lattice, and often only 
in the classical limit. A variational approach employing 
the Gibbs-Bogoliubov inequality fnjjllj can approximate 
the thermodynamics in all dimensions, either quantum or 
classical. Define a function $ which bounds the exact free 
energy function F(T) from above: 

F(T) < $ (a t , T)=F + (H- H ) , (2) 

where Tlo is a trial Hamiltonian which depends on ad- 
justable parameters {a,i} used to minimize the right hand 
side of the expression (||). F is the free energy of 
the trial Hamiltonian. The average <>o is taken with 
respect to Ho- If H.Q is harmonic, then minimization 
leads to temperature-dependent frequencies of TIq. This 
procedure is also known as the "self-consistent phonon 
method." We choose the trial Hamiltonian 




All quantities are dimensionless; ye is a displaced posi- 
tion variable xg + (— 1) u. If one raises T, the dimeriza- 
tion amplitude u goes to zero at a transition temperature 
T c ; H.q (Eq. |^) describes uncoupled oscillators with fre- 
quency w 2 , = Wq(1 + 4i5 2 cos 2 (fcet/2)). We study here only 
the classical limit for comparison with the MD simulation 
and the TM solution. The values of loq, 5 and u which 
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minimize (g) depend on £ = K3/K1 and on the inverse 
temperature /3. Let $mf denote the minimum value of 
the trial free energy $. From this we can compute the 
specific heat in MFA, C M f = -Td 2 <S> MF /dT 2 . We define 
a displacement correlation function G as 



G(m) 



l£<(-l)V-l)*^*fm) 



(4) 



The square of the order parameter is the correlation func- 
tion at large distance G (m — > 00). In MFA it is the vari- 
ational parameter u 2 . The correlation length, defined as 
decay rate of the correlation function at large m, is given 
in MFA by 
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Using the phonon dispersion law of the trial Hamilto- 
nian (g) we define a density of phonon states X>mf(w) = 
(a/2ir)dk/duj, which is a function of temperature- 
dependent parameters loq and S. 

Another quantity of the interest is the function 



F(x) 




x t + xt+t) 



(6) 



which gives the distribution of the separation be- 
tween nearest atoms. At the zero temperature xi = 
u(— and distribution is two delta functions F(x) = 
(S(x — 2u) + S(x + 2u)) /2. At non-zero temperature it 
can be evaluated in MFA 



F M f(x) 
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III. TRANSFER MATRIX METHOD. 

The TM method has been applied to one dimensional 
anharmonic systems |l^,[l4|]. To evaluate the classical 
partition function, one first writes it as 

N 



Z = jj J dzi exp {-(3f {z i+ x,Zi)) (8) 

i— 1 

/ (z i+ i, Zi) = - \z 2 +1 + \zf +1 + | (zi+i - Zif (9) 



where Zi is (— l) l Xi and obeys periodic boundary condi- 
tions zjv+i = z\. Integration of Z x S (z! — z±) over the 
variable z' simplifies the process when the delta function 
is expanded in a complete set of functions 



6(z'- Zl ) =5>* CO *„(*!) (10) 

n 

where the functions ^ n {z) satisfy the integral equation 
dz' exp (-/?/ (z, z')) *„ (z') = exp (-/3e„) * n (z) (11) 



Then the answer for the partition function in the thermo- 
dynamic limit TV — > 00 is Z = exp (~A/3e t ), where St is 
the "ground state" (exp (— /3e t ) is the largest eigenvalue) 
of the integral equation ([ll]). The leading term of the 
displacement correlation function is |L2] 



G(m) 
L 



{9't\y\9t) | 2 exp(-mZ c ) 
1 



(12) 
(13) 



where (and e' t ) are the first excited state of the integral 
equation ([ll]) . Then l c is the correlation length. 

We make the kernel of the integral equation symmetric 
by transforming the wavefunction: 



T(z) = *(z)exp 



--z z + -z- 
2 4 



(14) 



To find the smallest e„ we use a variational trial function: 



T tr (z) = exp 



4 , y 2 2 

— z H z 

4 2 



3=0 



(15) 



Here we have introduced a set of adjustable parameters 
61, 62, {ckj} used to minimize e% T (61, 62, {ay}) defined as: 



exp (-/3e t r) 



/ / dzdz'T tI .( 2 )K(z,z')T tr (2') 



/^T? r (z) 

X(z, z') = exp (/(z, z') + /(z', z))^) (16) 

In the decoupled case £ = the trial wave function (|l5| ) 
is exact with 61 =62 = (3/2 and m = 0. To find the first 
excited state we use an orthogonal function T( r = a;T tr . 
Then we minimize e' tr which is defined similarly to e t r in 
expression (|l6|), but with T tr replaced by T' tI . First we do 
minimization of (^|) only with respect to the parameters 
bi and b 2 setting all a.j to zero except ao = 1. Then we 
fix the values of 61, 62 and minimize equation (^) with 
respect to the set of {<x,-}, j = l..m. By increasing the 
number of adjustable parameters m we find the converged 
value of £t = min{e t r}. Usually m = 10 is sufficient to 
get an accurate answer with a relative difference less than 
10~ 10 between subsequent steps. 

Knowing the free energy e t we can deduce the specific 
heat by taking its second derivative with respect to tem- 
perature. Using the fact that Et is an extremum with 
respect to the adjustable parameters, we can evaluate its 
first derivative by computing the partial derivative with 
respect to temperature of the right hand side of Eq. (|l^) . 
Then a second derivative is found from finite differences 
of the computed first derivatives. 
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IV. NUMERICAL SIMULATION. 

We solved Newton's equations for the chain of atoms 
governed by Hamiltonian ([j]) . The system was a chain of 
250 atoms with periodic boundary conditions. A small 
time step dt = 0.014 was used. Each simulation lasted 
for time 70,000 (5 x I0 e dt). The position and veloc- 
ity of each atom were updated using the Verlet algo- 
rithm jl3] with error proportional to dt 4 . To control 
the temperature, the atoms suffered random collisions 
with "gas molecules" of equal mass. A collision was 
modeled by exchanging an atom's velocity with a gas 
molecule's velocity, which was randomly generated ac- 
cording to the Maxwell-Boltzmann distribution. Colli- 
sions occurred randomly but on an average of once every 
10 dt; the atom that suffered the collision was also chosen 
randomly. Atoms were initialized to be in the positions 
(xi = (— l) £ uo), with random initial velocities. The spe- 
cific heat was calculated using 

C= <E*>-<E>> (17) 

Each simulation was divided into 4 periods of 17,500 time 
units (1,250,000 dt). Over each of these 4 periods, the 
total energy E in the system was measured at 100,000 
randomly chosen time intervals. The specific heat was 
calculated for each of the 4 periods. The standard devi- 
ation of the 4 calculations served as an error bar. The 
correlation function G(to), for m = to 30, was calcu- 
lated at 400,000 randomly chosen time intervals and then 
averaged. The value of — ln(G(m)) was graphed versus 
to. The range of to values appropriate for linear fitting 
was chosen by eye. The slope of the linear fit was taken 
as the reciprocal correlation length. 

The frequency spectrum was similarly calculated from 
the velocity correlation function 



1 N 



(18) 



It was necessary to calculate from t = up to 3000 
in order for x(i) to approach 0. Each value of x(t) was 
calculated 1249 times and then averaged. At the end of 
simulation, the frequency spectrum was calculated by a 
Fourier transform 



x(t) cos(ujt)dt 



(19) 



The nearest neighbor distribution function (0) was found 
by averaging histograms for the displacement distribu- 
tion over 1000 randomly chosen times. 



V. RESULTS AND DISCUSSION. 

The exact solution has no ordered state at any T > 0. 
In MFA there is a critical temperature T c at which the 



trial free energies evaluated at u ^ and u = coincide, 
indicating a first order transition. The resulting order 
parameter u 2 versus temperature is shown on the top 
graph of Fig.[j]. The qualitative behavior is similar for all 
£, so we show results only for £ = 0.5. 

The specific heat in MFA is shown on the bottom graph 
of Fig|l] along with TM and MD solutions. The numerical 
MD and TM specific heats agree well. The uncertainty in 
the MD specific heat arises from statistical error, which 
should diminish as N^ 1 ^ 2 , when we take an average over 
a large number N of time slices. 
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FIG. 1. Specific heat, inverse correlation length and order 
parameter versus temperature for the case £ = Kz/ki = 0.5. 
The solid lines are exact solutions by the TM method; 
long-dashed line is MFA; diamonds are MD results. TM and 
MD results agree. 

There is a big discrepancy in the exact and MF solu- 
tions for specific heat, except in high and low T limits, 
where they coincide. To understand why MFA is accu- 
rate in the extreme T limits, consider the nearest neigh- 
bor distribution function in MFA and MD cases plotted 
on Fig. H for different T. At low T atoms are not far 
from the bottom of the well xg — (— 1) Uq. In this situa- 
tion approximation of the actual Hamiltonian ([!]) by the 
harmonic version (|3|) works well. In the opposite limit of 
high T harmonic approximation also works well, because 
the coupling term dominates the total energy. Averages 
of the first and the second terms of the MF Hamiltonian 
(§), namely u 2 < yt > 2 /2 and uj 2 S 2 < y t + y e+1 > 2 /2, 
coincide at 6 2 — 3/4, which happens at T = 0.74. At 
T = 0.8, exact and MFA results agree well for both spe- 
cific heat and nearest neighbor distribution function. Fig. 
U shows that the maximum deviation in F{x) happens at 
T = 0.28, just above T c = 0.26. The numerical result 
tells us that particles are mostly in the potential wells, 
while the MFA has maximum probability when particles 
are in undisplaced central sites. 

Results for the MF, TM, and MD correlation lengths 
are plotted on Fig.|l|. The latter two coincide. The error 
bars of the MD solution come from the least square fit 
of the numerical correlation function. These error bars 
are smaller than the size of the points and not shown 
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on Fig|l]. The correlation length diverges at T = and 
fits the ID Ising solution Aexp (B/T), where coefficients 
(A, B) depend on £, and approximately equal (0.18, 0.6) 
for £ = 0.5. The MFA gives a correlation length smaller 
than the exact one at all T and gives wrong behavior 
at T < T c . Above the MF T c the predicted dl c /dT is 
consistent with the exact solution. 



to the double well barrier Tb = 0.25. Particles with en- 
ergy > Tq take a very long time to overcome the barrier, 
which results in a contribution to the density of states at 
low frequencies. When we raise T well above T = 0.25, 
this low frequency contribution disappears, and the spec- 
trum becomes more similar to the MFA result. 
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FIG. 2. Nearest neighbor distribution function in MFA and 
MD solution for T=0.08 (solid line), T=0.16 (dashed line), 
T=0.28 (dotted line), T=0.8 (long-dashed line). Atomic dis- 
placement is measured in units of (ki/k.2) 1 ^ 2 ■ 

The maximum specific heat, which depends on £, oc- 
curs in MFA at T c . The temperature T max of the specific 
heat maximum in the exact solution is w 50% smaller 
than the T c of MFA, but tracks the MF transition point. 
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FIG. 3. Density of the phonon states versus frequency 
(in units of (ki/M) 1/2 ) in MFA and MD for T=0.08 (solid 
line), T=0.16 (dashed line), T=0.28 (dotted line), and T=0.8 
(long-dashed line). 

The density of phonon states found by MD and by 
MFA are shown on Fig. ||. In the low T limit again 
we see agreement between the MFA and MD solutions. 
When T is raised, the phonon density in MFA becomes 
broader and moves to lower frequency. Since the trial 
Hamiltonian (0) docs not have translational symmetry, 
no acoustic modes occur; beyond a lower bound, the spec- 
trum has no states. According to MD, a non-zero density 
of states appears at low frequency when T is comparable 
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